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Abstract 

The need for accurate material models to simulate the deformation, damage and failure of polymer 
matrix composites under impact conditions is becoming critical as these materials are gaining increased 
usage in the aerospace and automotive industries. While there are several composite material models 
currently available within commercial transient dynamic finite element codes, several features have been 
identified as being lacking in the currently available material models that could substantially enhance the 
predictive capability of the impact simulations. A specific desired feature pertains to the incorporation of 
both plasticity and damage within the material model. Another desired feature relates to using 
experimentally based tabulated stress-strain input to define the evolution of plasticity and damage as 
opposed to specifying discrete input properties (such as modulus and strength) and employing analytical 
functions to track the response of the material. To begin to address these needs, a combined plasticity and 
damage model suitable for use with both solid and shell elements is being developed for implementation 
within the commercial code LS-DYNA. The plasticity model is based on extending the Tsai-Wu 
composite failure model into a strain-hardening based orthotropic plasticity model with a non-associative 
flow rule. The evolution of the yield surface is determined based on tabulated stress-strain curves in the 
various normal and shear directions and is tracked using the effective plastic strain. The effective plastic 
strain is computed by using the non-associative flow rule in combination with appropriate numerical 
methods. To compute the evolution of damage, a strain equivalent semi-coupled formulation is used, in 
which a load in one direction results in a stiffness reduction in multiple coordinate directions. A specific 
laminated composite is examined to demonstrate the process of characterizing and analyzing the response 
of a composite using the developed model. 
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Introduction 


As composite materials are gaining increased use in aircraft components where impact resistance 
under high energy impact conditions is important (such as the turbine engine fan case), the need for 
accurate material models to simulate the deformation, damage and failure response of polymer matrix 
composites under impact conditions is becoming more critical. Within commercially available transient 
dynamic finite element codes such as LS-DYNA (Ref. 1), there are several material models currently 
available for application to the analysis of composites. The available models include relatively simple 
models such as a derivative of the Chang-Chang model (Ref. 2), where criteria related to ratios of stresses 
to failure strengths are used to signify failure, and the composite elastic constants are selectively reduced 
based on the failure mode. More sophisticated material models incorporated within transient dynamic 
finite element codes include continuum damage mechanics based models such as the model developed by 
Matzenmiller, et al. (Ref. 3). In this approach, the initiation and accumulation of damage is assumed to be 
the primary driver of any nonlinearity in the composite response. The failure stresses and strains of the 
material in each of the coordinate directions are specified by the user, the evolution of the damage is 
computed based on “damage variables”, and the nonlinearity of the material stress-strain response is 
approximated based on the input data and the evolution of the damage. While not necessarily included 
within commercial codes currently, another approach that has been used to model the strain rate 
dependent response of a composite is to assume that all of the nonlinearity is due to deformation 
mechanisms. An example of this approach was developed by Sun and Chen (Ref. 4), where a quadratic 
plastic potential function was developed and the plastic strains were computed based on the gradient of 
the plastic potential function. Stress-strain curves obtained for various fiber orientation angles are used to 
characterize the coefficients in the plastic potential function based on the values required to collapse the 
various separate stress-strain curves into a master curve. 

While the material models discussed above and other models have been utilized with some level of 
success in modeling the nonlinear and impact response of polymer composites, there are some areas 
where the predictive capability can be improved. Most importantly, the existing models often require 
significant a priori knowledge of the damage and failure in the analyzed structure such that their use as 
predictive tools can be limited. While these models generally assume that the nonlinear response of the 
composite is due to either deformation mechanisms (such as plasticity) or damage mechanisms, an 
improved model should have the capability to simulate the actual material behavior in which the material 
nonlinearity is due to a combination of both deformation and damage mechanisms. The input to current 
material models generally consists of point-wise properties (such as a specified failure stress or failure 
strain) that lead to curve fit approximations to the material stress-strain curves. This type of approach 
leads either to models with only a few parameters, which provide a crude approximation at best to the 
actual stress-strain curve, or to models with many parameters which require a large number of complex 
tests to characterize. An improved approach would be to use tabulated data from a well-defined set of 
experiments to accurately define the complete stress-strain response of the material. Furthermore, many 
of the existing models are only suitable for use with two-dimensional shell elements, which cannot 
capture the through-thickness response, which may be significant in impact applications. Ideally, a fully 
three-dimensional formulation suitable for use with solid elements would be desirable, along with a shell 
element formulation. 

To begin to address these needs, a multi-institution consortium has been formed to develop and 
implement a new composite material model within LS-DYNA, which will be implemented as MAT 213. 
The material model is meant to be a fully generalized model suitable for use with any composite 
architecture (unidirectional, laminated or textile). For the deformation model, the commonly used Tsai- 
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Wu composite failure criteria has been generalized and extended to a strain-hardening model with a 
quadratic yield function and a non-associative flow rule. The coefficients of the yield function for the new 
composite model are determined based on tabulated stress-strain curves in the various normal and shear 
directions, along with selected off-axis curves. The non-associative flow rule is used to compute the 
components of the plastic strain along with the effective plastic strain. The evolution of the yield stresses 
in the various coordinate directions is tracked based on the current value of the effective plastic strain. For 
the damage model, a strain equivalent formulation has been developed, which allows the plasticity and 
damage calculations to be uncoupled, and the plasticity calculations to take place in the effective 
(undamaged) stress space. In traditional damage mechanics models such as the one developed by 
Matzenmiller et al. (Ref. 3), a load in a particular coordinate direction is assumed to result in a stiffness 
reduction only in the direction of the applied load. However, for reasons to be discussed later in this 
paper, in the current model a semi-coupled formulation is developed in which a load in one direction 
results in a stiffness reduction in all of the coordinate directions. 

In the following sections of this paper, a summary of the derivation of the rate-independent 
deformation model is presented. The procedures to be used to characterize the material constants in the 
yield function and the flow law are discussed. Next, the rationale for, and detailed derivation of, the semi- 
coupled damage model is discussed along with a summary of the test matrix that will be required to 
properly characterize and validate the developed model. Finally, selected verification studies that have 
been conducted to ensure that the deformation model has been implemented correctly are presented. 


Deformation Model Derivation 

A general quadratic three-dimensional orthotropic yield function based on the Tsai-Wu failure model 
is specified as follows, where 1, 2, and 3 refer to the principal material directions. 
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In the yield function, ay represents the stresses and Fy and Ft are coefficients that vary based on the 
current values of the yield stresses in the various coordinate directions. By allowing the coefficients to 
vary, the yield surface evolution and hardening in each of the material directions can be precisely defined. 
The values of the normal and shear coefficients can be determined by simplifying the yield function for 
the case of unidirectional tensile and compressive loading in each of the coordinate directions along with 
shear tests in each of the shear directions, with results as shown below 
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In the above equation, the stresses are the current value of the yield stresses in the normal and shear 
directions (determined using procedures to be discussed below), where the superscript '/’indicates the 
tensile yield stress, and the superscript C indicates the absolute value of the compressive yield stress. To 
determine the values of the off-axis coefficients (which are required to capture the stress interaction 
effects), the results from 45° off-axis tests in the various coordinate directions can be used. An important 
point to note is that due to experimental or numerical variability, or alternatively just due to the 
fundamental behavior of the material, computing the off-diagonal terms of the yield function in this 
manner may result in a yield function that is not convex (which is a requirement for plasticity theory 
(Ref. 5)). As a result, to satisfy the requirements of the chosen yield function, the off-diagonal terms may 
need to be adjusted based on the values of the other coefficients in the yield function in order to ensure 
convexity of the yield surface. 

A non-associative flow rule is used to compute the evolution of the components of plastic strain. The 
plastic potential for the flow rule is shown below 


h = 

I ( 3 ) 
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where o, ; are the current values of the stresses and Hj are independent coefficients, which are assumed to 
remain constant. The values of the coefficients are computed based on average plastic Poisson’s ratios 
(Ref. 6). The plastic potential function in Equation (3) is used in a flow law to compute the components 
of the plastic strain rate, where the usual normality hypothesis from classical plasticity (Ref. 5) is assumed 

to apply and the variable X is a scalar plastic multiplier. 
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Given the flow law, the principal of the equivalence of plastic work (Ref. 5) can be used to determine 
expressions for the effective stress and effective plastic strain. By following this procedure, one can 
conclude that the plastic potential function h can be defined as the effective stress and the plastic 
multiplier can be defined as the effective plastic strain rate. 

To compute the current value of the yield stresses needed for the yield function, the common practice 
in plasticity constitutive equations is to use analytical functions to define the evolution of the stresses as a 
function of the components of plastic strain (or the effective plastic strain). Alternatively, in the 
developed model tabulated stress-strain curves are used to track the yield stress evolution. The user is 
required to input twelve stress versus plastic strain curves. Specifically, the required curves include 
uniaxial tension curves in each of the normal directions (1,2,3), uniaxial compression curves in each of 
the normal directions (1,2,3), shear stress-strain curves in each of the shear directions (1-2, 2-3 and 3-1), 
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and 45° off-axis tension curves in each of the 1-2, 2-3 and 3-1 planes. The 45° curves are required in 
order to properly capture the stress interaction effects. By utilizing tabulated stress-strain curves to track 
the evolution of the deformation response, the experimental stress-strain response of the material can be 
captured exactly without any curve fit approximations. The required stress-strain data can be obtained 
either from actual experimental test results or by appropriate numerical experiments utilizing stand-alone 
codes. Currently, only static test data is considered. Future efforts will involve adding strain rate and 
temperature dependent effects to the computations. To track the evolution of the deformation response 
along each of the stress-strain curves, the effective plastic strain is chosen to be the tracking parameter. 
Using a numerical procedure based on the radial return method (Ref. 5) in combination with an iterative 
approach, the effective plastic strain is computed for each time/load step. The stresses for each of the 
tabulated input curves corresponding to the current value of the effective plastic strain are then used to 
compute the yield function coefficients. 

Damage Model Derivation 

The deformation portion of the material model provides the majority of the capability of the model to 
simulate the nonlinear stress-strain response of the composite. However, in order to capture the nonlinear 
unloading and local softening of the stress-strain response often observed in composites (Ref. 7), a 
complementary damage law is required. In the damage law formulation, strain equivalence is assumed, in 
which for every time step the total, elastic and plastic strains in the actual and effective (undamaged) 
stress spaces are the same (Ref. 8). The utilization of strain equivalence permits the plasticity and damage 
calculations to be uncoupled, as all of the plasticity computations can take place in the effective stress 
space. 

The first step in the development of the damage model is to relate the actual stresses to a set of 
effective stresses by use of a damage tensor M 


o - M a e ff (5) 

The effective stress rate tensor can be related to the total and plastic strain rate tensors by use of the 
standard elasto-plastic constitutive equation 


® eff ~ ^p) (^) 

where C is the standard elastic stiffness matrix and the actual total and plastic strain rate tensors are used 
due to the strain equivalence assumption. By differentiating Equation (5) and substituting in Equation (6), 
the actual stress rate can be written in terms of the total and plastic strain rates and the actual stress. In the 
following expression, Voigt notation is assumed to be appropriate and the damage tensor is assumed to be 
invertible. 


o = M<j e ff + M G e ff 


(7) 


An algorithm to carry out the uncoupled plasticity/damage analysis is summarized below. In the 
algorithm, the superscript represents values computed in the previous time step, and the superscript 
“m+ 1” indicates values to be computed in the current time step. In the first step of the algorithm, the actual 
stresses are converted into effective stresses using the damage tensor M. In the second step, the plasticity 
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calculations are carried out in the effective stress space to compute the current value of the plastic strain 
rate, and the effective stress values are updated. Next, in step 3 the damage tensor is modified based on 
the computed plastic strain rate. Finally, in step 4 the modified damage tensor is used to compute the 
updated values of the actual stresses based on the updated effective stresses. The algorithm is summarized 
symbolically below, where At is the time step. 

1. o« // = (m ->)V 

2 - (8) 

3. M' ,+1 =M"+Am[e /) ] 

4. a n+l = 

To ensure that the strain equivalence assumption is valid for this approach, the equations utilized in the 
algorithm should be able to be manipulated to yield an equation similar to that shown in the second 
expression in Equation (7), which was derived based solely on the strain equivalence assumption. By 
starting with the equation shown in step 4 of the algorithm, substituting the equation used in step 2 of the 
algorithm in for the modified effective stress, and breaking up the modified damage tensor into the sum of 
the current value of the damage tensor and the increment in the damage tensor (AM), the following set of 
equations is developed. 
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To a first order approximation, the last expression in Equation (9) corresponds to the second expression in 
Equation (7), which indicates that strain equivalence assumption is appropriate for the combined 
plasticity/damage algorithm. 

Given the usual assumption that the actual stress tensor and the effective stress tensor are symmetric, 
the actual stresses can be related to the effective stresses in the following manner, where the damage 
tensor M is assumed to have a maximum of 36 independent components. 


CGx ' 


( a eff^ 

V XX 

a yy 


a eff 

°yy 

® zz 
^ xy 

= [M] 

G eff 

U ZZ 

G e ff 

^xy 

Qyz 


c eff 

®XZ 7 


eff 

V a xz j 


( 10 ) 


N ASA/TM— 20 15-218849 


6 



In many damage mechanics models for composites, for example (Refs. 3 and 7 ), the damage tensor is 
assumed to be diagonal or manipulated to be a diagonal tensor, leading to the following form. 

M n 0 0 0 0 0 

0 M22 0 0 0 0 

0 0 AT33 0 0 0 

000 m 44 0 0 

0 0 0 0 M 55 0 

0 0 0 0 0 m 66 

The shear terms M44, M55 and Mm, can be independent (such as in Matzenmiller, et al. (Ref. 3 )) or 
functions of the normal damage terms Mn, M21 and M3 3. 

The implication of a diagonal damage tensor is that loading the composite in a particular coordinate 
direction only leads to a stiffness reduction in the direction of the load due to the formation of matrix 
cracks perpendicular to the direction of the load. However, several recent experimental studies 
(Refs. 9 , 10 , and 11 ) have shown that in actual composites, particularly those with complex fiber 
architectures, a load in one coordinate direction can lead to stiffness reductions in multiple coordinate 
directions. One example of this phenomena can be seen in Figure 1 , which came from research conducted 
by Salem and Wilmoth (Ref. 11 ). The figure shows a triaxially braided composite which has been loaded 
in tension in the transverse direction and then reloaded in compression in the longitudinal direction. As 
can be seen in the figure, after loading in the transverse direction significant damage above and beyond 
simple matrix cracks was present in the material. As a result, when the material was reloaded in the 
longitudinal direction, the measured longitudinal modulus was significantly reduced from the baseline 
value, indicating that the damage resulting from the transverse load affected the stiffness in the 
longitudinal direction. There have been limited attempts to incorporate this damage coupling into an 
analytical technique, for example by Voyiadjis and Park (Ref. 12 ) and Bednarcyk, et al. (Ref. 13 ). 
However, these efforts were developed within the overall concept of a damage mechanics theory in which 
all of the nonlinearity of the composite response was assumed to be due to damage mechanisms. 
Furthermore, in the developed theoretical approaches analytical functions were used to track the evolution 
of the damage and the reduction of the stiffness. The efforts described in the current paper are geared 
towards developing a damage theory, uncoupled from the plasticity theory, which tracks the stiffness 
reduction and damage accumulation as a function of the plastic strain by the use of tabulated input. 

One approach to incorporating the coupling of damage modes would be to use a non-diagonal damage 
tensor, such as the one shown below for the case of plane stress. 

'O \Mu M n M 16 faf" 

Gyy = M 2 \ M 2 i M 26 ( 12 ) 

k® y _~^61 Mq 2 M 56 - ^ rt\T ^ 

However, while this formulation would allow for directional coupling, it would have the side effect of a 
unidirectional load in the actual stress space resulting in a multiaxial load in the effective stress space. For 
the strain equivalent combined plasticity-damage formulation envisioned for this model, this would be an 
undesirable side effect as the plasticity calculations could be adversely affected due to the introduction of 
nonphysical stresses. 
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Longitudinal Load 



Figure 1. — Damage observed in triaxially braided composite 
subjected to transverse loading followed by longitudinal 
loading. 


I dxx // dxx 

[ xx h xx 

^ XX 



Figure 2. — Coupled damage resulting from loading in x coordinate direction. 

To avoid the undesired stress coupling, a diagonal damage tensor is required. However, to account for 
the damage interaction in at least a semi-coupled sense, each term in the diagonal damage matrix should 
be a function of the plastic strains in each of the normal and shear coordinate directions, as follows for the 
example of the M\ \ term for the plane stress case 

M\\ =Afn(e&,e£,,e£,) (13) 

To explain this concept graphically, a schematic is shown in Figure 2 for the case of loading in the x 
coordinate direction. A plastic strain is applied to an undamaged specimen, with an original area A xx 
perpendicular to the x axis and an original area ^ perpendicular to the y axis. The undamaged modulus 
in the x direction is E xx and the undamaged modulus in the y direction is equal to Eyy. The specimen is 
damaged due to the plastic strain. The original specimen is unloaded and reloaded elastically in the x 
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direction. Due to the damage, the reloaded specimen has a reduced area in the x direction of A^ x and a 
reduced modulus in the x direction of E $£ x . The reduced area and modulus are a function of the damage 
induced by the plastic strain in the x direction as follows 


pdxx — (i _ rfxxlpP r 
^ xx ~ y u xx\'xx)y-‘xx 

Adxx — 1 _ r/xxlp P M 
™XX — V {A ‘XX\'XX)) rL XX 


where is the damage in the x direction due to a load in the x direction. Alternatively, if the damaged 
specimen was reloaded elastically in the y direction, due to the assumed damage coupling the reloaded 
specimen would have a reduced area in the y direction of Ay f, x and a reduced modulus in the y direction 

of Eyy X due to the load in the x direction. The reduced area and modulus are again a function of the 
damage induced by the plastic strain in the x direction as follows 


^^-dyy^jEyy 

Ap^(\-dyy(zP x )]A yy 


where d is the damage in the y direction due to a load in the x direction. Similar arguments can be 
made and equations developed for the situation where the original specimen is loaded plastically in the y 
direction. 

The next issue is how to properly model the damage coupling for the case of multiaxial loading. A 
simple case for consideration is a composite specimen being simultaneously strained plastically in the x 
and y directions. One way to approach the problem would be to assume that the loss in area in the x 
direction due to the straining in the x direction is as follows. 

A — Adxx — A x/xx n fit 

Act Act — AoAxt V 10 / 

Next, one can assume that the loss in area in the x direction due to the loading in the y direction is as 
follows 


A _ A^yy — A A xx 
Act Act ~ ^xx u yy 


( 17 ) 


where A^ y is the reduced area in the x direction resulting from a load in the y direction and d xx is the 

damage in the x direction resulting from a load in the y direction. The total loss in area in the x direction 
could be computed by adding together the two separate area losses to obtain a total loss in area. 


2A XX 




( 18 ) 


The ratio of the total damaged area to the undamaged area in the x direction can now be computed, which 
theoretically would lead to the total amount of damage in the x direction in the composite under biaxial 
loading 
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( 19 ) 


A X x-2A xx + 
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The relationship between the actual stress and effective stress in the x direction would then be the 
following. 


= (l -d xx is p )-d xx iz p 

\ A u xx \pxx ) u yy \*yy JJ KJ xx 


( 20 ) 


The error in this approach is that the assumption is made that the load in the x and y directions are 
both acting on undamaged areas. In reality, the loads are acting on damaged areas, and just adding to the 
damaged area. For example, if one loaded the material in the y direction first, the reduced area in the x 

direction would be equal to A x p and the reduced modulus in the x direction would be equal to E x p . If 
one would then subsequently load the material in the x direction, the baseline area in the x direction 
would not equal the original area A xx , but the reduced area A x p . Likewise, the baseline modulus in the x 

direction would not be equal to the original modulus E xx , but instead the reduced modulus E x g . 
Therefore, the loading in the x direction would result in the following further reduction in the area and 
modulus in the x direction. 

Eip 4- dp(pp))E*p’ 4- dg[ 4J (i - </4 e 4kvx 

A$P = (l- = (l - dg{ z p xx )) (l - 

These results suggest that the relation between the actual stress and the effective stress should be based on 
a multiplicative combination of the damage terms as opposed to an additive combination of the damage 
terms. For example, for the case of plane stress, the relation between the actual and effective stresses 
could be expressed as follows 


= (l - dg )(l - dg )(l - dg )ag 

Oyy = \-dS%-d»%-d%\4 ( 22 ) 

where for each of the damage terms the subscript indicates the direction of the load which initiates the 
particular increment of damage and the superscript indicates the direction in which the damage takes 
place. Note that for the full three-dimensional case the stress in a particular coordinate direction is a 
function of the damage due to loading in all of the coordinate directions (x, y, z, xy, xz and yz). By using 
a polynomial to describe the damage, the coupled terms represent the reduction to the degree of damage 
resulting from the fact that in a multiaxial loading case the area reductions are combined. 

To properly characterize the damage model, an extensive set of test data is required. Due to the 

tabulated nature of the input, each of the damage parameters ( dp , dp , etc.) has to be determined as a 

function of the plastic strain in a particular coordinate direction (such as eg ). For example, to determine 
the damage terms for the case of loading in the x direction, a composite specimen has to be loaded to a 
certain plastic strain level in the x direction. The material is then unloaded to a state of zero stress, and 
then reloaded elastically in each of the coordinate directions to determine the reduced modulus of the 
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material in each of the coordinate directions. Expressions such as those in Equations (14) and (15) can 
then be used to determine the required damage parameters for the particular value of plastic strain. The 
process needs to be repeated for multiple values of plastic strain in the x direction in order to establish the 
full characterization of the variation of the damage parameters as a function of the plastic strain in the x 
direction. 


Verification Studies For Deformation Model 

A set of verification studies for the deformation portion of the material model were conducted using 
data for a T800S/3900-2B unidirectional composite (Ref. 14). This particular composite is composed of 
intermediate modulus, high strength fibers embedded within a toughened epoxy matrix. Full details of the 
verification studies will be provided in a future paper, however a summary is provided here for 
completeness. An important point to note is that at the current time only the deformation portion of the 
material model has been implemented numerically within the LS-DYNA computer code, so only 
deformation analyses will be discussed here. Future efforts will involve implementing and verifying the 
damage theory described above. Flowever, since the evolution of the parameters in the damage model will 
be based on the plastic strains, ensuring that the nonlinear deformation response of the composite is being 
properly simulated by the material model is critical to ensuring that the subsequent characterization and 
verification of the damage model is correct. 

The input data utilized for the verification studies is a combination of actual experimental data obtained 
by Raju and Acosta (Ref. 14) and numerical simulations. The numerical simulations were conducted to 
obtain required stress-strain curves which were not available from the provided experimental data. To 
conduct the numerical simulations, micromechanics analyses were conducted where the properties of the 
fiber and matrix were used to simulate the overall response of the composite. The micromechanics analyses 
were conducted by using a combination of purely analytical simulations conducted using the NASA Glenn 
developed MAC/GMC code based on the Generalized Method of Cells (Ref. 15) and finite element models 
where the fiber and matrix were explicitly simulated. The fiber was assumed be linear elastic with 
transversely isotropic properties and the matrix was assumed to be isotropic with an elastic -plastic material 
response. A fiber volume fraction of 0.54 was assumed based on data presented in Bogert, et al (Ref. 16). 
The constituent material properties were chosen to correlate with composite experimental data obtained by 
Raju and Acosta (Ref. 14) and Bogert, et al. (Ref. 16). The elastic modulus and Poisson’s ratio of the matrix 
were chosen to be values representative of epoxy based materials (Ref. 17). The fiber longitudinal modulus, 
transverse modulus and longitudinal Poisson’s ratio were chosen such that the composite longitudinal 
modulus, transverse modulus and longitudinal Poisson’s ratio computed by the micromechanics analyses 
correlated with the values obtained by Bogert, et al. (Ref. 16). For the fiber transverse Poisson’s ratio, a 
value representative of that utilized for similar carbon fibers (Ref. 1 7) was used. The fiber in-plane shear 
modulus was chosen such that the computed composite in-plane shear modulus correlated to the composite 
in-plane shear modulus determined by Raju and Acosta (Ref. 14). The transverse shear moduli were 
computed by appropriate assumptions based on the assumed transverse isotropy of the unidirectional 
material. The yield stress for the matrix was selected such that the computed in-plane shear stress-strain 
curve for the composite roughly correlated with the experimental stress-strain curve determined by Raju and 
Acosta (Ref. 14). Full details of the numerical experiments and the process used to correlate the properties 
for the numerical experiments will be provided in a future paper. The correlated fiber and matrix properties 
used for the analyses are provided in Table 1. 

To conduct the verification studies, finite element models such as the ones shown in Figure 3 for 
tension and Figure 4 for shear were used. Note that 64 eight noded solid elements were used for the 
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analyses. For the tension simulations, the nodes on the left hand side of the model were constrained and a 
constant displacement in the x direction was applied to the nodes on the right face. For the shear 
simulations, the nodes on the bottom surface were constrained and a displacement in the x direction was 
applied to the top surface. 


TABLE 1.— FIBER AND MATRIX CONSTITUTIVE 
PROPERTIES USED FOR MICROMECHANICS ANALYSES 


Property 

Fiber 

Matrix 

£u,GPa 

275.6 

3.45 

E 22 , GPa 

15.5 

15.5 

Vl2 

0.20 

0.35 

V23 

0.25 

0.35 

G 12 , GPa 

103.4 

1.27 

CJyield, MPcl 

N/A 

137.8 




Figure 4. — Finite element model for shear verification analyses. 
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Simulated stress-strain curves were computed for a variety of load cases, examples of which are 
shown in Figure 5 for the case of [0°] tension, Figure 6 for the case of in-plane shear in the x-y plane and 
Figure 7 for [45°] off-axis tension in the x-y plane. In all of the figures the “Experimental” curves were 
the curves that were provided as input data for the material model which were based on the actual or 
numerically generated stress-strain curves, and the “Simulated” curves were the curves computed using 
the material model. In all cases, the simulated curves correlated very well to the input curves, indicating 
that the material model was able to accurately represent the nonlinear deformation of the composite in 
terms of correctly replicating the input data. Future efforts will involve conducting more complex 
validation analyses which will include analyzing the deformation response of more complex laminated 
composites as well as simulating the deformation response of a composite under impact conditions. 



0.000 0.005 0.010 0.015 0.020 

Strain 

Figure 5. — Verification analyses for 0° tension test. 
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Stress (MPa) Stress (MPa) 


1-2/3-1 Shear Test 



0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 


Strain 

Figure 6. — Verification analyses for in-plane shear test. 
Off-Axis Tension (1-2 plane, 45 Degree) 



0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050 

Strain 

Figure 7. — Verification Analyses for 45° off-axis tension test. 
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Conclusions 


A generalized composite model suitable for use in polymer composite impact simulations has been 
developed. The theory for the rate independent deformation and damage portions of the composite model 
have been developed and numerical implementation of the deformation model has been completed. The 
complete composite model will be implemented into the LS-DYNA commercial transient dynamic finite 
code as MAT 213. For the deformation model, the Tsai-Wu composite failure model has been generalized 
into an orthotropic yield function with a non-associative flow rule. Tabulated stress-strain data is utilized 
to track the deformation response of the material, using the effective plastic strain as the tracking variable. 
A strain equivalent damage model has been developed in which loading the material in a particular 
coordinate direction can lead to damage in multiple coordinate directions. The actual and effective 
stresses are related by multiplicative combinations of the various damage variables. 

Future efforts will involve generalizing the deformation model to incoiporate the ability to simulate 
the effects of strain rate and temperature on the material response. The damage model will be numerically 
implemented within the LS-DYNA code. Methods to model failure and element removal will be 
developed and implemented into LS-DYNA. Extensive additional verification and validation studies will 
be conducted to verify the accuracy and capability of the overall material model. Overall, when 
completed the material model when implemented into MAT 213 will provide significant improvements to 
the state of the art in the modeling of the impact response of polymer composites. 
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